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Abstract: We present the finite temperature spectra of both bottomonium and char- 
monium, obtained from a consistent lattice QCD based potential picture. Starting 
point is the complex in-medium potential extracted on full QCD lattices with dynam¬ 
ical u,d and s quarks, generated by the HotQCD collaboration. Using the generalized 
Gauss law approach, vetted in a previous study on quenched QCD, we fit Re[U] with a 
single temperature dependent parameter the Debye screening mass, and confirm 
the up to now tentative values of Im[U]. The obtained analytic expression for the 
complex potential allows us to compute quarkonium spectral functions by solving an 
appropriate Schrodinger equation. These spectra exhibit thermal widths, which are free 
from the resolution artifacts that plague direct reconstructions from Euclidean corre¬ 
lators using Bayesian methods. In the present adiabatic setting, we find clear evidence 
for sequential melting and derive melting temperatures for the different bound states. 
Quarkonium is gradually weakened by both screening (Re[U]) and scattering (Im[U]) 
effects that in combination lead to a shift of their in-medium spectral features to smaller 
frequencies, contrary to the mass gain of elementary particles at finite temperature. 
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1 Introduction 

Heavy quarkonium, the bound states of a heavy quark and anti-quark (cc or 66), are 
a unique tool to simplify the complexity inherent in the physics of the strong inter¬ 
actions. Instead of having to deploy the full force of quantum field theory, we may 
consider a Schrodinger equation with an effective potential to describe the dynamics of 
these bound states in vacuum and in-medium. At zero temperature the two main char¬ 
acteristics of QCD, asymptotic freedom and conhnement manifest themselves clearly 
in the form of a Coulombic behavior with a running coupling at small distances and 
a non-perturbative linear rise at large distances [1]. he. we can learn about key fea¬ 
tures of the strong interactions by inspecting this simple potential. Vice versa, from 
the knowledge of the potential, we can reproduce even quantitatively a significant part 
of vacuum quarkonium physics (e.g. the bound state spectrum below the open heavy 
quark threshold [2, 3]). 
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The question of how such a potential arises from the microscopic theory of QCD has 
been answered in detail by advances made in the research on the effective held theories 
(EFT) NRQCD and pNRQCD. The latter has been introduced at hrst at T = 0 in 
[4] and investigated in the non-perturbative regime in [5, 6 ], with both cases being 
reviewed in detail in [7]. The generalization to hnite temperature in a perturbative 
setting followed in [ 8 ]. 

In the presence of a separation of scales, the combination of NRQCD and pNRQCD 
offers a systematic power counting prescription in the heavy quark velocity n, to setup a 
simplihed description of the bound state evolution at soft energy scales (Egoft ~ rngv) 
in terms of singlet and color octet wavefunctions. he. we do not have to explicitly 
describe the physics at the much higher hard scale (Ehard ~ 2 mQ)h And indeed the 
heavy quark rest mass (in this work we use rric ~ 1.472GeV and = 4.882GeV) and 
the characteristic scale of quantum huctuations in QGD, usually denoted by Aqcd — 
200MeV are far enough apart to warrant a quantitatively reliable potential description 
in vacuum. 

Effective held theory has furthermore reminded us that the potential between two 
inhnitely heavy quarks is an inherently Minkowski-time quantity. Starting out from 
fundamental QGD, the evolution of a QQ can be described by the thermally averaged 
unequal-time point-split meson-meson correlator 

= <^M(x,y,t)M^(xo,yo,0)^, (1.1) 


with the meson operator dehned as M(x,y,t) = Q(x,f) 7 '^f/[x,y]Q(y,f). 7 ^ refers to 
the Dirac matrices and 7/[x, y] denotes a straight Wilson line connecting (x, f) and 
(y, t). The relative coordinate enters as r = x — y. 

In the limit of inhnite quark mass, where the static constituents are truly test 
particles, ( 1 . 1 ) can be identihed with the medium averaged unequal time correlation 
function of quarkonium singlet wavefunctions in the language of the EFT. It is this 
wavefunction correlator for which a Schroedinger equation and thus an in-medium 
potential will be established. For static quarks the spatial separation r = |r| becomes 
an external parameter of the theory and their evolution traces out a rectangle over 
time. Formally it has been shown that eq.(l.l) reduces to the rectangular real-time 
Wilson loop [7] 


Wn{r,t) 


Tr (exp 


-ig 



dx>^Ay 

H' 



( 1 . 2 ) 


^The strength of the EFT approach lies in the fact that any residual influence of the physics of the 
hard scale can be systematically included via the matching of Wilson coefficients, which goes beyond 
the ability of direct methods, such as e.g. the historic T — 0 Wilson loop approach [9] 
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As such, this object still contains the physics of all scales: hard, soft and ultra-soft, i.e. 
in general it fulhlls an equation of motion [10] 

idtWu{r,t) = ^{r,t)Wu{r,t) (1.3) 


with <h(r, t) being a time- and space dependent complex function. 

In case that a potential picture for the static QQ system is applicable, the function 
<I>(r, t) has to asymptote to a time independent function V{r). In turn it will dominate 
the evolution of Ihn(r, t) at times much larger than the intrinsic scales of e.g. the gluons 
and light quarks in the medium. Formally we may then write 


V{r) 


zdtW{t,r) 

hm ——-— 

t^oo W{t,r) 


(1.4) 


This limit reflects the fact that in order to replace a retarded, i.e. gluon mediated 
interaction with an instantaneous potential, the gluons must have been exchanged 
between the heavy quarks multiple times. 

In the presence of the additional scale T we need to ascertain how reliable a de¬ 
scription of realistic, i.e. hnite-mass, quarkonium is in terms of the static potential 
alone. Two forms of so called relativistic corrections can adversely affect the accuracy 
of the lowest order approximation. The hrst kind remains fully within the potential 
picture, i.e. hnite mass corrections to the static potential can become signihcant, as 
they are proportional to the relative velocity of the heavy quarkonium system. These 
corrections can be systematically computed (see e.g. [1, 5, 6]), in vacuum they are 
found to be small and we expect the same to be true at hnite temperature. Their 
non-perturbative determination from hnite temperature lattice QCD will be the aim of 
a future study. 

The other contributions are so called non-potential ehects, arising from physics that 
cannot be recast into the form of a time independent term in the Schrodinger equation. 
One way they can enter in the perturbative formulation of ehective held theories if the 
ultra-soft gluons are still kept as explicit degrees of freedom (see e.g. [7, 11]). In a non- 
perturbative setting, if non-potential ehects become sizable, the late time evolution of 
the Wilson loop cannot be described by a simple time-independent potential V{r). As 
will be discussed in the next section we hnd that at the temperatures and inter-quark 
distance investigated here, the size of such contributions remains insignihcant. 

Interestingly the dehnition of the potential from the Wilson loop in real-time coin¬ 
cides with the late r limit in Euclidean time in the case of T = 0. This constitutes the 
basis for the successful extraction of the static zero temperature potential from lattice 
QCD simulations, which are carried out solely in an imaginary time setting. At hnite 
temperature, the real-time dehnition of the potential does not change, however the 
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imaginary time axis becomes compactified and its finite extend encodes vital physics 
information, i.e. the inverse temperature. Hence the straight forward connection be¬ 
tween the late Minkowski time limit and the Euclidean Wilson loop at maximum t = f3 
is absent. 

For more than two decades, the absence of an effective-field theory based definition 
for the the in-medium potential has led theorists to embrace model potentials that 
were defined directly from Euclidean time observables [12], readily calculable in lattice 
QCD. In particular two quantities have gained popularity as model potentials, the color 
singlet free energies in Coulomb gauge 


g/3F(l)(r) 


Tr 





r2(r) = exp 


-ig r dTA;(:r,T)T‘ 


(1.5) 


defined from the correlator of Polyakov loops r2(r) and a derived quantity the color 
singlet internal energies — TS. While these quantities exhibit a behavior 

compatible with the expectations for e.g. Debye screening of the interaction between 
the heavy quarks in the deconfined phase, it could be shown that they do not match 
the potential for the quark anti-quark system at finite temperature [15, 16]. 

And indeed neither one of these fully real quantities by themselves can take the role 
of a static in-medium heavy quark potential, as we have learned from a ground breaking 
series of works starting with Laine et. ah [8, 10] in 2007. In their study the real-time 
definition (1.4) was evaluated in a resummed perturbative framework, called the hard- 
thermal loop approximation. This approach contains a gauge invariant resummation 
of an infinite number of Feynman diagrams and has been shown to capture many key 
features of QCD at high temperature reliably. As the Wilson loop in Minkowski time 
is an in general complex quantity, the authors observed that the potential too takes on 
complex values at finite temperature 


Chtl(?’) — —ds 




rriu + 


+ iT(l){m^r) 


+ 0(9“), 


( 1 . 6 ) 


with 


0(a;) = 2 



dz 


z 

(^2 + 1)2 



(1.7) 


Here we absorb a factor Cp m. the definition of the coupling constant to 

connect to the conventions in the phenomenology literature. It was furthermore shown 
that the real-part of this complex potential itself does not coincide with the color- 
single free energies in hard-thermal loop resummed perturbation theory in [15, 16], 
even though the absolute deviations are comparatively small. 
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The fact that the in-medium potential is a complex quantity not only reflects a 
quantitative change but necessitates a qualitatively different perspective on the physics 
it describes, he. besides screening of the force between the heavy quarks (Debye 
screening) seen in the real-part, the effects of scattering of light medium degrees [10, 17] 
with the heavy quarks (Landau damping) related to Im[l/] further weaken the bound 
state in the QCD medium. In fact, depending on the hierarchy of scales present, the 
imaginary part can be related to different phenomena, such as the breakup of a color 
singlet to an octet conhguration [8] and in turn to the dissociation of the QQ system 
into gluons [18]. 

Only these effects taken together can give a consistent picture of heavy quark 
bound states at hnite temperature, he. a simple estimate of the dissolution of a heavy 
quarkonium state solely on the basis of vanishing binding energy is not sufficient. This 
in turn has consequences for heavy quarkonium phenomenology in general, where e.g. 
the melting temperatures enter as input into transport model calculations. 

We would like to stress that the in-medium potential discussed here does not di¬ 
rectly govern the evolution of the actual wave-function of the heavy quarkonium system. 
By construction it instead enters in the Schrodinger equation of the real-time Wilson 
loop, which in the EFT language corresponds to the thermally averaged correlator of 
unequal time wavefunctions. The imaginary part of the potential hence describes the 
decay of this correlator over time, which does not directly imply the annihilation of 
the heavy-quarks. In fact, due to the non-relativistic approximation we operate under, 
Q and Q remain in the system forever. However, even if the norm of the quarkonium 
wavefunction is preserved in unitary time evolution, its correlation with the initial state 
can still decay with time, a phenomenon known as decoherence. It is an active area 
of research to answer the question how to connect the complex in-medium potential to 
the evolution of the bound state wave-function, which has not yet been answered in 
the effective held theory setting of pNRQCD. The concept of open-quantum systems 
(see e.g. [19]) has proven insightful in this regard. 

One possible way to elucidate the physics encoded in the complex potential is to 
compute the spectral function p{u;) of heavy quarkonium. This quantity is related 
to the heavy quark current-current correlator, which can be obtained from (1.1) by 
carefully taking the limit of vanishing point splitting. In section 3 we will hence solve 
an appropriate Schrodinger equation to compute p{u)). Once computed one can observe 
how the formerly delta-like bound state peaks present at T = 0 broaden and shift as 
screening and scattering modihes the QQ state with increasing temperature. Choosing 
a popular criterion of melting temperature [10, 60] at the point where binding energy 
and spectral width coincide we can furthermore determine the point of dissolution for 
different bottomonium and charmonium states. 
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Changes in spectral structure are of particular interest as they are directly related 
to changes in the dilepton emission from quarkonium decay. The dilepton emission rate 
is given by a simple product of the in-medium spectral function with the Bose-Einstein 
factor [13] 


dRu 

d^P 


37r3p2 


nB(Pii)p(P'), 


( 1 . 8 ) 


where Qg denotes the electric charge of the heavy quark considered in units of e, the 
four momentum is P = (po, p) and the hnite mass of the leptons has been neglected (for 
a more detailed discussion of the above formula see ref. [14]). he. if a bound state or its 
remnant appears as well dehned peaked feature in an in-medium spectral function, the 
area under such a structure informs us about the experimentally accessible dilepton 
emission of that state. Note that while the above relation is only applicable to the 
decay of a quarkonium state in a thermalized static plasma, it can nevertheless provide 
us with vital physics insight as we will see in sec. 4.1. 

This paper is organized as follows. We start sec. 2 with a review on recent progress 
in the extraction of the values of the complex in-medium potential from Euclidean time 
lattice QCD simulations. In order to deploy the lattice potential in phenomenological 
applications, sec. 2.2 discusses the generalized Gauss law ansatz developed in ref. [30], 
which provides an analytic formula for Re[V] and Im[V] depending on a single temper¬ 
ature dependent parameter mn the Debye mass. Tuning mz) we are able to reproduce 
the lattice values of Re[V] and conhrm the yet tentative values for Im[V]. These values 
for m£) are compared to perturbation theory. Section 3 is concerned with using the 
continuum corrected potential to investigate the phenomenology of in-medium heavy 
quarkonium. We will calculate the spectra for the bottomonium and charmonium S- 
wave channel and investigate their modihcation with increasing temperature. Besides 
giving the melting temperatures for different states, we will determine the T' to J/T ra¬ 
tio at the chiral crossover and give a rough estimate for the suppression of bottomonium 
in a heavy-ion collision compared to p -|- p. We close with a conclusion in sec. 5. 


2 Lattice QCD potential and Debye screening mass 

While the perturbative computation of the potential has contributed signihcantly to 
our understanding of the physics involved, it is not sufficient for the description of 
the experimentally relevant temperature regime around the phase transition. There 
the quark gluon plasma can indeed be considered strongly interacting, exemplihed, 
e.g. by the large value of the trace anomaly [20]. This calls for an evaluation of the 
potential dehnition (1.4) in lattice QCD, which at hrst seems unfeasible, since direct 
access to real-time quantities, such as the Wilson loop (1.2) is not possible. Conceptual 
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and technical advances in the extraction of real-time information from lattice QCD 
simulations over the last few years however have made such an evaluation possible, as 
we will discuss below. 

2.1 The complex in-medium potential from lattice QCD 

The real-time Wilson loop cannot be directly accessed in lattice QCD simulations, as 
they are performed in Euclidean time. One strategy, proposed in [21] and applied for 
the hrst time in [22] to circumvent this problem is to resort to a spectral decomposition 
of the Wilson loop, which simply amounts to a Fourier transform over a positive dehnite 
spectral function p 



/ 


due = Wu{t,r). 


due ‘^'^puioj^r) o 


The fact that the time dependence is explicit in the integral kernel tells us that both the 
real-time and Euclidean time Wilson loop are described by the same spectral function. 
Hence extracting the values of p from imaginary time simulation data will give access 
to the real-time Wilson loop and in turn to the potential. Inserted into the dehning 
equation (1.4) the spectrum itself can be related to the values of the potential 



( 2 . 1 ) 


Unfortunately extracting the spectrum from Euclidean time simulation data is an 
inherently ill-dehned inverse problem, as one seeks to determine the form of a continuous 
function from a hnite and noisy set of individual points. Only by using additional 
prior information, such as the positive dehniteness of the spectrum or smoothness 
assumptions is it possible to give meaning to this problem, a strategy usually referred 
to as Bayesian inference. Established implementation of this approach, such as the 
Maximum Entropy Method [23] or extended MEM [24] have been shown to fail to 
produce satisfactory results when deployed in the extraction of spectra from Wilson 
loops [25], and it needed the development of a new Bayesian reconstruction prescription 
[26] before quantitatively robust results were obtained. Even with this new method the 
spectrum can be determined only in a certain range of frequencies, given by the energies 
resolved on the lattice, which makes a brute force evaluation of (2.1) impossible. A 
selection of reconstructed spectra for the two extreme cases of the lowest and highest 
investigated temperatures around Tc are shown in Eg. 1. 

This further difficulty can also be overcome from a careful inspection of the spectral 
structure of the Wilson loop. Just as we did to arrive at eq. (1.4) let us assume that 
a potential picture is valid, i.e. the late time evolution of lTn(f, r) is dominated by 
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Figure 1. (top) A selection of spectral functions from the Euclidean Wilson Line correlators 
in Coulomb gauge, extracted from 48^ x 12 HotQCD asqtad lattices using a novel Bayesian 
reconstruction prescription [26] at the lowest (left) and highest (right) temperatures. To 
better observe their peak and shoulder structures, individual spectra are shifted manually in 
frequency. For T /Tc = 0.86 a well defined lowest lying peak is present at any of the shown 
distances r, while at T/Tq = 1.66 the reduced physical extend degrades the reconstruction 
at larger r. (bottom) Illustration of the skewed Lorentzian (SL) nature of the individual 
spectral peaks. Using eq.(2.2) up to quadratic polynomial order, we fit the points in the dark 
red region close to the maximum in each case. The resulting fit functions (dashed orange) 
reproduce the shape of the spectra (blue circles) quantitatively far beyond the actual fitting 
region. Note that this indicates that the peak is a Lorentzian and not e.g. a Gaussian 

a time independent function lim^^oo‘^’(^, ^) = V{r). As was proven in [27], in this 
case the Wilson loop spectrum will contain a well dehned lowest lying peak of skewed 
Lorentzian form embedded in a polynomial background 

lImU(r)lcos[Re(7oo(r)] - (ReU(r) - u;)sin[Recroo(r)] 

^ ImU(r)2 + (ReU(r) — a;)2 

+ Co(r) + Ci(r)(Rel/(r) — w) + C 2 (r)(Rel/(r) — _ 

The position and width of this peak are related to the real- and imaginary part of the 
potential respectively. The skewness and background contributions on the other hand 
are identihed with remnants of the physics at scales above the soft-scale. 
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We can now tnrn this relation aronnd and nse it as a non-pertnrbative criterion for 
the applicability of a potential description. If the Wilson loop spectrnm does contain 
a well dehned lowest lying peaked featnre of skewed Lorentzian type, it will lead to a 
late-time evolntion governed by a well dehned potential. As shown in the examples of 
reconstrncted spectra in £g. 1 we do indeed find snch well prononnced peaked featnres 
and trongh a £t establish that they of a skewed Lorentzian form^. This leads ns to 
conclnde that an in-medinm potential description is warranted and the valnes of V(r) 
can be extracted via the application of eq.(2.2). 

In practice, we look for the lowest lying peak in the spectra reconstructed from the 
lattice QCD observables, £t it around the full width at half maximum with the skewed 
Lorentzian form (2.2) and read off the value of Re[ld] and Im[I/] encoded in it. This 
extraction strategy laid out above has been successfully tested using hard-thermal loop 
perturbation theory [25], where both the Wilson loop in Euclidean time, its spectrum 
and the corresponding potential are known. 

In this study we use the values of the heavy quark potential [28] extracted^ from 
full QCD lattices generated by the HotQCD collaboration [29], containing dynamical 
u,d, and s quarks. The medium quarks on these lattices with spatial extend W = 48 
are described by the asqtad action and are tuned to lie on a line of constant physics 
with a physical strange quark mass and light u,d quark masses close to their physical 
values ms/mu,d = 20 (For more details see [29]). In addition to the values of the 
potential around the phase transition at T /Tq = 0.68 — 1.66, which had been extracted 
in a previous study [28], we add here the values from two additional low temperature 
ensembles close to T ~ 0, which will be used to calibrate the analytic expression for the 
temperature dependence of the potential in the next subsection, tab. 1 summarizes the 
relevant simulation parameters for our ensembles and gives the temperatures in relation 
to the chiral crossover transition temperature on these lattices at Tc = 172.5 MeV. 

Since the hnite temperature lattices only contain W = 12 lattice points in Eu¬ 
clidean time direction, we expect that while a robust determination of spectral peak 
positions is possible, the accuracy for spectral widths will be insufficient to make reli¬ 
able statements about Im[I/]. The values obtained for the real-part are given by the 
colored points in the left panel of £g. 2, the tentative values for Im[I/] from the spectral 

^If the peaks were instead Gaussian pg(^) = coexp [—(oj —m)^/r^], inserting them into the defining 
formula (2.1) would lead to a divergent expression, since / doj ujpG{oj)e~^‘^* / f doj pG{uj)e~^‘^* =m — 
shows an unphysical linear increase in Im[V] over time [22]. 

^Note that the observable we use here to extract the potential is not the Wilson loop but the Wilson 
line correlator in Coulomb gauge. The reason is that the latter does not contain cusp divergences that 
make evaluating the Wilson loop very costly on the lattice. We have checked that the result remains 
gauge invariant by relaxing the Coulomb condition and applying random gauge transformations. 
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a [fm] 

0.1 

0.057 

0.111 

0.1 

0.09 
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0.071 

0.068 

0.057 

48^ X Nr 

48 

48 

12 

12 

12 

12 

12 

12 

12 

N 

4 ’ meas 

350 

163 

1295 

1340 

1015 

1270 

1220 

1150 

1130 


Table 1. Parameters of the isotropic HotQCD 48^ x Nr lattices [29] with asqtad action 
{mi = ms/20,Tc ~ 172.5MeV) used in this study. 


widths are plotted as lightly colored points in the background of the right panel. 

We will postpone the discussion of the temperature dependence of the potential to 
the end of the next subsection. 

2.2 Gauss law parametrization of the potential 

In order to investigate heavy quarkonium spectra by solving a Schrodinger equation 
with the proper complex heavy-quark potential, we require an easily evaluable expres¬ 
sion for both Re[ld] and Im[l/] at a given temperature. To this end we have recently 
proposed a field theoretically motivated functional form, that depends on a single tem¬ 
perature dependent parameter, the Debye screening mass [30]. 

The idea behind this approach is as follows. The physics of heavy quarkonium at 
zero temperature can be described in a potential picture with two distinct characteris¬ 
tics, a Coulombic part at small distances and a linearly rising part at large distances^. 
Hence if we wish to understand the in-medium modihcation of the inter-quark poten¬ 
tial we need to understand how a test charge associated with either one of these two 
distinct potentials will react to the presence of medium charge carriers in its surround¬ 
ings. To describe the effect of these charges we will use the hard-thermal loop medium 
permittivity, in essence making the ansatz of a test particle with a particular electric 
field configuration being immersed in a bath of weakly interacting quarks and gluons. 

The basis for the derivation of the analytic expression for the complex potential 
lies in the generalized Gauss law derived in [31] 

=47rg5(r-0. (2.3) 

It is applicable to a point charge with (color) electric held E = which covers 

both the Coulombic a = —l,q = cxs, [ds] = 1, as well as string-like potential a = 1, g = 

^Note that by allowing the linear term to contribute down to the shortest distances we partially 
absorb the running of the strong coupling 
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a, [a] = GeV^ relevant for heavy quarkonium. In case of a Conlombic potential we 
can incorporate the in-medinm effects in a straight forward manner by transforming 
(2.3) into Fonrier space and mnltiplying the right hand side with a (possibly complex) 
permittivity, taken here from hard-thermal loop pertnrbation theory 


e \p,mD) 


P 


+ rn^ 


— zttT 


pw?jy 


(p2 -|- m 


ir-‘ 


(2.4) 


Transforming back to coordinate space yields 


— V^V'c(r) +mj^Vc{r) = as{^Ti5{f) — iTmjjip^mDr)^, 


(2.5) 


with 


p{x) = 2 / dp 


sin(pa;) p 


px p‘ 


+ 1 ’ 


( 2 , 6 ) 


which constitntes an integro-differential eqnation with linear-response character for 
the in-medinm modihed Conlomb potential. The strength of the medinm effects is 
clearly governed by the parameter mD, which we have proposed as a non-pertnrbative 
dehnition of the Debye mass in ref. [30]. Solving (2.5) with the appropriate bonndary 
conditions ReVc(’")lr=oo ~ ~ ^ = 0 reprodnces the 

pertnrbative potential of Laine et. al [10] given in eq.(1.6). 

In the case of a string-like potential, transforming eq.(2.3) into Fonrier space is 
impractical. Instead we rely on the assnmption that the in-medinm charge distribntion 
fonnd in the Coulomb case is the same as for the linearly rising potential i.e. a property 
of the medium and not of the test charge. As was shown in [30] the strength of the 
medium effects in the resulting linear-response type equation is governed not by the 
Debye mass alone but instead by the parameter /i^ = mjj-?- 


1 / \ 

- ^ = cr(^47r(5(r) - iTm]jp{mDr)j. (2.7) 

Solving for the real part of the medium modihed string potential leads to an analytic 
expression in terms of parabolic cylinder functions Dy{x) 

K^V.{r) = -^^D_,{V2^,r) + ^^ ( 2 . 8 ) 

P ^ [4J P 

The imaginary part on the other hand can be written in a closed form using the 
Wronskian method, which yields 

ImVs(r) = -i - ^'ip{pr) = -ia^Ttl^^pr), (2.9) 
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Figure 2. (left) Real-part of the in-medium heavy quark potential on Nf = 2 -|- 1 asqtad 
lattices (points) with fits that establish the values of the Debye mass (solid line). Errorbands 
denote changes from varying the value oi rriD within its fit uncertainty, (right) prediction of 
the imaginary part of the potential (solid curves), together with the tentative values (light 
points) extracted from the asqtad lattices with Nr = 12. 


with 


'ip{x) = D_i/ 2 {\/ 2 x) / dyReD_i/ 2 {iV 2 y)y‘^ip{ymD/n) 

Jo 

roo 

+ReD_i/ 2 {iV 2 x) / dy D_i/ 2 {V 2 y)y‘^Lp{ymD/y) 

J X 

poo 

-£>- 1 / 2 ( 0 ) / dy D_i/2{V2y)y'^ip{vmD/f‘)- 


( 2 . 10 ) 


Note that even though we have used a weak coupling ansatz for the permittivity, 
which is only appropriate at high temperatures, its insertion into on the Gauss law 
has produced an expression for the potential with linear-response character. In turn 
it seems to smoothly connect to zero temperature, as the value of the Debye mass 
can in principle be set to zero. In particular the real part of the Gauss law derived 
in-medium potential reduces to the T = 0 Gornell-potential at m^) = 0. This bodes 
well for applying the derived expression to htting the lattice QGD extracted potential 
values even at temperatures below the formal range of applicability of weak-coupling 
methods. 
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2.3 Determination of m^) from the lattice potential 

For a consistent determination of the single temperatnre dependent parameter rriD, we 
assnme that neither the strong conpling nor the string tension depend on the medinm 
temperatnre, i.e. all modihcation emerges from the snrronnding light qnarks and gln- 
ons. Thns we will need to £x the valnes of as, cr, as well as the arbitrary scale dependent 
constant shift of Re[ld] at T = 0. Since lattice QCD simnlations are performed on fi¬ 
nite size lattices at a hnite lattice spacing, they necessarily operate at a low bnt hnite 
temperatnre. In onr case we will hence nse the newly added data from the two en¬ 
sembles close to zero temperatnre at /3i = 6.9 and (32 = 7.48. We hnd that the valnes 
of the vacnnm parameters vary slightly between the two ensembles. Hence a linear 
interpolation is nsed at intermediate lattice spacings.^ As can be seen from the dark 
and light bine curves at the top of the right panel of £g. 2 the lattice values for Re[17] 
at low temperature indeed show a well pronounced Coulombic and linear behavior that 
is excellently reproduced by the fit parameters given in tab. 2. I.e. neither the running 
of the coupling at small distances nor logarithmic corrections at large distances are 
signihcant for the regime investigated here. 



/3=6.9 

/3=7.48 

Otg 

0.456 ±0.027 

0.385 ±0.006 

^ [GeV] 

0.470 ±0.01 

0.515 ±0.003 

c [GeV] 

1.760 ±0.034 

2.648 ±0.009 


Table 2. Values for the vacuum potential parameters from the low temperature lattice QCD 
ensembles 

With the vacuum values set and since the explicit expression for both Re[R] and 
Im[l/] derived above only depends on a single parameter, we continue by determining 
rriD from a £t to the real part of the lattice QCD extracted values alone. The resulting 
curves are given as solid lines in the right panel of fig. 2, the corresponding Debye 
masses are collected in tab. 3 and plotted in fig. 3. 

We hnd that tuning mo allows us to indeed reproduce both the qualitative and 
quantitative behavior of Re[R] without problem, if the interpolated T = 0 constants 
are used. The strength of the generalized Gauss law Ansatz is that it now allows us to 
predict the values of the imaginary part of the potential, which was deemed unreliable 
in a previous study due to the fact that on the hnite temperature lattices the Wilson 
lines were available at only twelve points in Euclidean direction. We are conhdent in 

®Note that the differences between the values of a at different (3 values might be related to an 
insufficiently precise setting of the scale for the asqtad lattices in [29] 
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6.8 

6.9 

7 

7.125 

7.25 

7.3 

7.48 

T/Tc 

0.86 

0.95 

1.06 

1.19 

1.34 

1.41 

1.66 

mo 

0.01(3) 

0.25(8) 

0.39(8) 

0.53(21) 

0.96(5) 

0.99(13) 

1.27(8) 

rriD 

T 

0.04(10) 

0.72(22) 

1.03(22) 

1.28(49) 

2.07(11) 

2.05(27) 

2.29(14) 


Table 3. Debye masses extracted from the isotropic HotQCD 48^ x 12 lattices with asqtad 
action. For use in phenomenology, a continuum corrected mo may be obtained from the ratio 
mj)/cr{(]) shown here, through a multiplication with the continuum value of a. 



Figure 3. (left, blue points) The normalized Debye mass {mo/y/<y) from a generalized Gauss- 
law fit to the real-part of the in-medium heavy quark potential on asqtad lattices, which we 
propose as input for phenomenological studies, (right), Blue points: temperature dependence 
of the Debye mass (m^/T) and in red a NLO HTL based fit of rriD 


the predictive capabilities of the approach, as it has been able to successfully reproduce 
Ini[l/] in a similar study in quenched QCD [30]. 

When plotting Im[l/], as implied by m£)(T) via eq. (2.5) and eq.(2.9) as solid lines 
in the right panel of £g. 2 the tentative data (light colored points) shows a surprisingly 
good agreement with these values at small separation distances r < 0.75fm down to 
T = 0.95T(7. Only at T = O.SOT^ we hnd large differences between the extracted 
and predicted values, which is probably due to the impossibility to resolve the tiny 
width of the spectral peaks at low temperatures using only = 12 data points. 
Since heavy quarkonium phenomenology in the quark gluon plasma requires knowledge 
about the potential only until the freezeout boundary is reached, i.e. slightly below the 
deconhnement transition, the apparent range of applicability of the £t seems to suffice. 

2.4 Continuum correction for the potential and mj) 

Several preparations are still in order, since for a meaningful phenomenological investi¬ 
gation, we need to use the continuum values for all parameters entering the potential. 
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In the absence of a true continuum extrapolation of our lattice results we resort to the 
following strategy. 

We £t the vacuum potential parameters, as, a and the constant term c by com¬ 
paring the energy levels of the corresponding Cornell-potential Hamiltonian to the 
experimentally measured masses of the bottomonium system, where we expect that 
hnite mass correction are insignihcant. More importantly for bottom quarks there 
exists a well controlled procedure to dehne their physical mass. A mass often de¬ 
ployed in pNRQCD effective theory calculations is the pole mass [32-34] (for bottom 
^poie ^ 4.93GeV) different from the usually quoted MS mass = 4.18). Its 

larger value reflects the fact that the quark mass should be evaluated at the soft scale, 
relevant for the physics of the bound state and not at the hard scale. A renormalization 
group flow towards lower energies is hence required, leading to the larger . More 
precisely, because of the so called renormalon ambiguity in the perturbative evaluation 
of the Wilson coefficients of the effective theory, one is lead to use the renormalon 
subtracted scheme [35]. he. one subtracts the renormalon ambiguity in the pole mass 
definition and reshuffles it into the constant part of the potential, which suffers from 
the same ambiguity. In this way, the ambiguities eventually cancel and we are lead 
to a quantitatively robust definition of the potential. In general also here one obtains 
larger values than those resulting from the MS scheme, such as the established 

= 4.882 ± 0.041 GeV (2.11) 

from ref. [35], which has been deployed e.g. in [7]. The errors are estimated to be as 
large as the highest calculated contribution of order C>(a|). 

This choice of mass allows us to reproduce the PDG values for the four S-wave states 
T(IS')—T(4S'), the averaged massed of the three known P-wave triplets Xhi)-P)—Xh{‘iP) 
as well as the lowest D-wave state T(lZl) to at least three significant digits if the 
remaining vacuum parameters are set to the following values 

c=-0.1767±0.0210GeV, d, = 0.5043±0.0298, = 0.415T0.015 GeV. (2.12) 

These continuum values are very close to the ones used in conventional quarkonium 
spectroscopy [2, 3] and are also compatible with our lattice data. Note that the effects 
of string breaking were introduced by hand, flattening off the potential at Vgb = 1.25fm 
[36]. A naive complex scaling analysis confirms that the four lowest S-wave states of 
the modified Gornell potential Hamiltonian indeed lie below the continuum threshold. 

Lattice QGD simulations are carried out in a finite box with a finite lattice spacing, 
hence discretization effects have to be taken into account. The former is related to the 
fact that the quarks do not take on their physical masses exactly and thus the chiral 
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crossover temperature on our lattices = 172.5MeV) is higher than the continuum 
value of T™"* = 155 ± 9MeV. In addition the change in lattice spacing leads to slightly 
different values of the vacuum parameters of the potential in our case, as discussed 
in sec. 2.3. Hence both effects should be corrected for. Ultimately, we will use the 
constants of eq.(2.12) together with the continuum corrected Debye mass 

= T/Tg^^^) = (2.13) 


where 


mp 


is given in tab. 3 and 


in (2.12). 


Let us take a closer look at the continuum corrected Debye masses and how they 
compare to perturbative estimates. In the work of ref. [37] it has been established that 
the Debye mass can be computed perturbatively only up to the leading order together 
with the logarithmic correction at next to leading order (NLO). The presence of a 
magnetic sector in QCD leads to the appearance of truly non-perturbative contributions 
to m^) at NLO, parametrized in the following by the two terms containing the constants 
ki,K 2 , which need to be determined from numerical simulations 


Wc I iU 

3 6 


^ , JN, ^ N, _ NJ'gipf -, 


+HiTg{pf + 


(2.14) 


To compare the temperature dependence of the Debye mass, we £x the non-perturbative 
constants Ki, K 2 while keeping /i = 2'kT constant. For the running of the coupling g{^) 
we utilize the four loop result of ref. [38] setting Aqcd = 0.2145GeV, appropriate for 
starting the renormalization group flow from a scale where Nf = 5 flavors are active. 

In the left panel of £g. 3 we show the values oi mo/T (blue points) together with 
the £t according to eq. (2.14). Even though the perturbative part of the formula 
for mD is applicable, if at all, at the highest temperatures investigated here, we End 
that the fit manages to pass through all values oi ttid even around the phase bound¬ 
ary. The non-perturbative contributions are non-negligible with ki = 0.84 ± 0.10 and 
^2 = —0.40 ± 0.03. Nevertheless the perturbation theory motivated £t (by chance) 
reproduces mr) down to temperatures, slightly below Tc and may hence be used to de¬ 
fine a phenomenological, lattice QCD validated, temperature dependence of the Debye 
screening mass. 


3 Quarkonium spectra 

In sec. 2 we managed to capture the functional form of both the real- and imaginary- 
part of the static potential by fitting a single temperature dependent parameter 
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the Debye mass, to Re[D] extracted on the lattice. We now take the next step and 
compute from it the in-medium spectral functions of the vector channel bottomonium 
and charmonium S-wave states at hnite mass by solving an appropriate Schrodinger 
equation. These spectra correspond to resting states with absolute momenta p = 0, 
similar to those usually investigated in direct lattice QCD or lattice NRQCD studies. 

The static heavy quark potential is a universal quantity, in the sense that it denotes 
the lowest order contribution in the non-relativistic expansion for both bottomonium 
and charmonium physics. Therefore we expect that the vacuum parameters htted in the 
bottomonium case do remain the same for the lighter flavor. An important difference 
between the two cases is that the binding energy of the charmonium ground state 
is close to Aqcd- This makes it impossible to set up a perturbative renormalization 
scheme for the charm mass, similar to the one we used for Bottom. Hence instead of 
calculating the renormalon subtracted mass, we use the charm mass as £t parameter 
and tune it to reproduce the masses of the stable S-wave states (J/T,T') and the 
averaged two P-wave triplets Xc{.^P)yXc{‘^P)- The resulting best fit value reads 

m^DGfit = 1.472 GeV. (3.1) 

Since hnite mass effects, in particular radiative corrections can become relevant for 
charmonium, we expect the agreement between the static potential based masses and 
the experimental T = 0 spectrum to be worse than for bottomonium. Indeed for the 
S-wave states and the IP triplet only agreement up to the second digit is found when 
these higher order effects are neglected. 

The vacuum parameters determined in sec. 2.4 and (3.1) constitute the basis from 
which we embark on the hnite temperature study, where all medium ehects on the static 
potential are summarized in the temperature dependence of the Debye mass parameter 
determined in sec. 2.3. We use the continuum string tension of the bottomonium ht to 
convert the values of mj:)/s/a of tab. 3, i.e. the right panel of hg. 3 for use in the contin¬ 
uum Schrodinger equation. Note that we have set up the potential parametrization in 
sec. 2 such that changing the Debye mass does not ahect the overall constant in Re[R]. 
I.e. at small enough separation distances, where temperature ehects are irrelevant, the 
values of Re[R] all agree independently of m^), which is expected of a correctly renor¬ 
malized potential and a necessary requirement for a meaningful interpretation of the 
quarkonium bound state physics at hnite T. tab. 4 and tab. 5 summarize the vacuum 
properties of several bottomonium and charmonium states that arise from the T = 0 
potential parameters. 
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states 

T(l^) 

T(2^) 

T(35) 

T(4^) 

XbilP) 

Xb{2P) 

Xbi^P) 

m [GeV] 

9.4603 

10.020 

10.353 

10.597 

9.92597 

10.269 

10.538 

[GeV] 

9.4603 

10.023 

10.355 

10.579 

9.88814 

10.252 

10.534 

(r) [GeV-i] 

1.489 

2.985 

4.385 

10.17 

2.435 

3.898 

5.586 

(r) [fm] 

0.2934 

0.5881 

0.8639 

2.004 

0.4797 

0.7679 

1.100 

m[GeV] 

1.1 

0.539 

0.206 

-0.038 

0.633 

0.29 

0.02 


Table 4. The masses, mean radii and distances to the BB threshold for bottomonium S-wave 
and P-wave states at T=0. 


states 

J/^(15) 

^'(2^) 

XcilP) 

Xc{2P) 

m [GeV] 

3.0969 

3.6717 

3.5089 

3.7918 

[GeV] 

3.0969 

3.6861 

3.4939 

3.9228 

(r) [GeV-i] 

2.861 

5.839 

4.136 

25.42 

(r) [fm] 

0.5635 

1.150 

0.814824 

5.00813 

m[GeV] 

0.639 

0.064 

0.227 

-0.056 


Table 5. The masses, mean radii and distances to the DD threshold for charmonium S-wave 
and P-wave states at T=0. 


3.1 Spectral functions from the Schrodinger equation 

All ingredients have been assembled for computing the vector channel spectral functions 
from the time evolution of the corresponding correlation function D^{t, r, r') governed 
by the complex in-medium potential®. In the following we deploy the Fourier space 
method developed in ref. [39], which solves 


H + i\lmV{T)\ 


D>(t,r,r') = idtD>(t 
D^(t, r, r') = idtD^(t 


with 


H 


2mQ - 


V" 

2mQ 


+ Re[I/](r) + 


and the starting condition 


D>(0,r,r') = - 


, r,r'), t>0 

(3.2) 

, r,r'), t<0 

(3.3) 

/(/ + 1) 

(3.4) 

mqr^ 


r'). 

(3.5) 


®Note again that it is not the Schrodinger equation for the quarkonium wavefunction but for the 
forward correlator that we are solving here. 
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For the correlator in frequency space we Fourier transform 


/ 


'OO 


D{uj,r,r') 


r, r'), 


(3.6) 


— OO 


from which the vector channel spectrum is obtained by taking the limit 



(3.7) 


Taking the point splitting to zero requires some effort in the practical implementation 
as discussed in appendix A of ref. [39]. 

In £g. 4 we give overview plots of the computed in-medium S-wave spectral func¬ 
tions for bottomonium and charmonium at several temperatures around the transition 
temperature. Note that the widths seen here are actual physical widths, related to 
hnite temperature effects, i.e. all bound states reduce to delta peak structures in the 
T = 0 limit. By construction, the T = 0 peaks coincide with the experimental values 
for the bottomonium states up to 3 digits and £t the two charmonium states below 
threshold up to few percents 

3.2 Properties of the in-medium spectral functions 

As expected from the difference in constituent quark mass, the bottomonium states 
react less severely to the medium at a given temperature compared to charmonium. 
In general we can observe in this adiabatic setting that the very narrow peaks at low 
temperature are affected in a hierarchical manner, the highest lying states, which are 
most loosely bound, begin to broaden and shift hrst, followed sequentially by the lower 
lying states. The combination of the effects of screening Re[V] and scattering Im[V] 
encoded in the potential lead to characteristic common changes for both quark flavors. 
Not only do the bound states broaden but also their mass shifts to lower values. This 
behavior is found not only for the ground state but for all the higher lying states before 
they eventually dissolve and become part of the continuum. 

To be more quantitative, if a narrow resonance pole lies close to the real frequency 
axis its spectrum can be described by a simple Breit-Wigner (BW). On the other hand 
for states that are close to melting, it is necessary to disentangle the remnant bound 
state signal from the continuum background. For such broad features a skewed Breit- 
Wigner needs to be considered [50], which reads 



^For T = 0, we actually add a small imaginary part in the potential to avoid having exact delta- 
functions in the spectrum, that would be impossible to visualize graphically. 
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Figure 4. (top panel) S-wave bottomonium in-medium spectral function based on a proper 
complex potential as defined in sec. 2.4. In this adiabatic scenario we find clear indications for 
sequential melting. T(4S') is almost gone already at 0.94Tc, T(3S') disappears rapidly close 
to Tc and T(25) starts to be washed out at 1.19Tc. Only T(IS') remains discernible at our 
highest lattice temperature l.bGTc*. (bottom panel) The charmonium in-medium spectral 
function based on the same complex potential as for bottomonium. In the absence of a 
well defined renormalon subtracted charm mass, its value is fitted to reproduce the known 
charmonium bound states at T = 0 giving rric = 1.472GeV. Also charmonium shows a 
sequential melting pattern, as we find that 'I'(25) disappear around Tq while J/T close 
to 1.41Tc. 


where E denotes the energy of the resonance, F its width and 5 the phase shift. 

Using the interpolated form for the Debye mass (2.14) we perform a temperature 
scan of the spectrum at every 6t = ^ MeV and £t the different peaks with the BW 
of eq.(3.8). From these hts we obtain the temperature dependence of the bound state 
width and mass of different quarkonium states, which are plotted in fig. 5 and fig. 6 
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Figure 5. Mass of the charmonium (left) and bottomonium (right) bound states from the 
position of their in-medium spectral peaks. The peak position decreases monotonously with 
temperature until the bound state disappears. The continuum threshold energy Fcont = 
Re[V(r —)• oo)] is shown as gray line. A hierarchial destabilization of the bound states in 
accord with the reduction in binding energy is seen. Note that depending on the skewness of 
the peak, the mass it encodes does not have to coincide with it’s apex position. The error 
bands reflect the uncertainty of the Debye mass determination. 

respectively. Another quantity particularly relevant to phenomenology is the area under 
each of the peaks, which we dehne here as the integrated area of the Breit-Wigner £t 
function A = and which is related to the total dilepton emission rate via eq.(1.8). 
Its values are given in £g. 7. We hnd that as temperature increases, the peak area 
at hrst remains rather constant, even if the peak becomes wider but eventually and 
abruptly begins to decrease rapidly towards zero. 

We can put these observations in the context of the in-medium modification of the 
potential shown in £g. 2. While the small distance part of the correctly renormalized 
Re[V] is virtually temperature independent, a signihcant screening of the linear rise 
at large distances occurs with increasing T. It is a peculiarity of the conhnement 
mechanism that with the diminishing remnant of the linear rise; also the threshold 
to the continuum i^cont is lowered. This is reflected in a decrease of the value of 
A'cont = Re[V(r —>■ oo)] (see the gray curve in fig. 5). In turn, the binding energy, 
defined from the difference between bound state mass and the continuum threshold 
energy is monotonously lowered. That is to say, the thermal fluctuations of the medium 
destabilize the bound state. Interestingly, once the threshold moves into the vicinity of 
a formerly hrmly bound quarkonium peak it pushes the spectral feature towards lower 
frequencies until it eventually disappears. 

Let us consider the mass shift observed in fig. 5. At first it might seem counterintu¬ 
itive, as the expectation for an elementary particle in a medium is exactly the opposite. 
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T[MeV] 


Figure 6. Width of the Charmonium (left) and Bottomonium (right) bound states, which 
increases monotonously with temperature. The error bands reflect the uncertainty of the 
Debye mass determination. 



T[MeV] 



Figure 7. Area under the bound-state peaks in the charmonium (left) and bottomonium 
(right) spectrum. Note the characteristic plateau region before a rapid decrease to zero sets 
in. The error bands reflect the uncertainty of the Debye mass determination. 


it will receive a thermal mass. For instance at LO, a single fermion receives a mass 
correction [40-42] of 

„2 2 ^ 2 ^ 

mQ{T) =mQ + 6mQ = mQ + ^^j^^, (3.9) 

and one might wonder if such a contribution should be added to our potential. As 
can be seen from (3.9) it is of higher order in the I/tuq expansion and hence does 
not contribute to the static potential. At the next order in the l/mq expansion, the 
potential V{r) is corrected by a term of the form Vi{r)/mQ [5]. If the two color charges 
are widely separated one expects that each of them obtains a thermal mass shift, hence 
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linir^oo On the other hand, for r —)■ 0 one can imagine the neutral 

bound state will not interact with the plasma so that Vi{r) = 0. The thermal mass shift 
(3.9) may hence be considered an upper bound on the thermal mass shift of the bound 
state. Since it is of higher order it is negligible in our non-relativistic computation, 
which can be checked using our value of as = Q^Cpj^'K and considering a temperature 
of T ~ 200MeV. One would obtain Sttiq = 7MeV for charm and Sttiq = 2MeV for 
bottom, insignihcant in comparison to the thermal shift observed in £g. 4, which easily 
reaches 100 — 300MeV (see also fig. 5). 

These hndings, based on a hrst principles lattice QCD based complex in-medium 
potential, are qualitatively similar to what had been observed in potential modeling 
studies that solve a Schrodinger equation with a potential with imaginary part inserted 
by hand (see e.g. [43]). There only the Coulombic contribution to Im[l/] was used, 
which leads to more stable behavior than in our case. The presence of the string- 
like vacuum potential contributes an additional term to Im[l/], which is of comparable 
size as the Coulombic contribution to Im[l/] at intermediate temperatures. On the 
other hand our results differ signihcantly from the spectra obtained in the T-matrix 
study of ref. [44]. There the authors use purely real model potentials which lead to 
in-medium states that appear to possess a significantly larger energy than their vacuum 
counterparts. Updated computations in that framework [45] are expected in the near 
future. Similar shifts of the resonance peaks to lower energies were also observed in a 
sum-rule based approach [46-48], but the magnitude of the shift is somewhat weaker 
there. A more thorough comparison of our results to direct studies in lattice QCD is 
part of ongoing work. 

For T ~ 250MeV one can also compare to the perturbative estimates of [49] and 
qualitatively good agreement between the spectra is found. The lattice spectra show 
slightly narrower peaks than the perturbative ones, due to the linearly rising part of the 
potential, which increase the binding energy. This effect is not captured in perturbation 
theory. On the other hand, the string effects vanishes quickly at high temperature and 
in addition the imaginary part of the lattice potential is larger due to string effects. 
Combined it seems to compensate the stronger binding through an increase in the width 
of the state. One central beneht of the lattice computation is the possibility to connect 
the high and low temperature spectra, which is cannot be realized in perturbation 
theory. 

A qualitative difference between the perturbative results and the lattice ones is the 
way the peak positions change as function of the temperature. In perturbation theory, 
which contains essentially only the Coulomb term, the bound state peaks move slightly 
to higher values of frequency [39] as T increases whereas we see here a clear decreases 
of the peak frequencies. 
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Figure 8. Charmonium (left) and Bottomonium (right) spectrum at T = l.lDTc with 
manually changed values for the imaginary part of the potential (0.1 ImV (black), ImV (red) 
and 2 ImV (blue)). While the peak position hardly changes unless the peak is close to melting, 
the width significantly depends on the strength of ImV). 


As the determination of the imaginary part of the potential represents a significant 
sonrce of nncertainty in this work, we close this section by stndying its effect on the 
spectrnm in more detail. To do so, we vary the strength of the imaginary part mnlti- 
plying it by a r-independent factor see £g. 8. The main effect of the imaginary part is 
to broaden the peak, withont changing its position and area, unless the peak is already 
close to the continuum, as is e.g. the case with T(3S') in fig. 8. That said, it is obvious 
that the corresponding states do melt more quickly with a large Im[l/], as the width 
reaches the binding energy much more quickly. 


4 In-medium quarkonium penomenology 

4.1 Charmonium at freezout 

At current heavy-ion colliders, such as RHIC and LHC, with collision energies above 
^/saa > 200GeV, the success of the statistical model of hadronization, to predict the 
yields of heavy quarkonium, supports the idea that charmonium completely dissolves 
in the created plasma. In turn essentially all charmonium bound states we observe in 
experiment would be generated via recombination that takes place at the freeze-out 
boundary, usually located slightly below the crossover temperature. In such a setting 
the ratio between the yields for J/T and T(2S') can be estimated from the difference 
in area under the corresponding peak structures in the in-medium spectral functions 
in a straight forward way. 
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w[GeV] w[GeV] 

Figure 9. charmonium (left) and bottomonium (right) spectra at Tc- To visualize the 
uncertainties of the computation three curves according to the values oi mo and ± 6m d 
are shown. 


Let us assume that due to the dynamical nature of the expanding hreball the freeze- 
out happens close to the chiral crossover temperature Tc- When inspecting the spectra 
at this temperature, as done in fig. 9 we find that there indeed exist two peaked structure 
related to the in-medium J/\k and fl''. he. we can compute the amount of in-medium 
dimuons produced by each of the two states via the product of the spectral function 
and Bose-Einstein distribution, integrated over the frequency range corresponding to 
each of the states via eq. (1.8). 

This however is not what is measured in experiment, since the charmonium states 
do not decay into leptons within the plasma but given their life-time, instead long after 
the plasma is diluted away. We should thus in principle project the states corresponding 
to the peaks of the hnite T spectra onto the T = 0 states. As no agreed upon method 
exists to do so, we here assume that the states inside the peaked structures will become 
real J/T or T' after freeze out and subsequently decay outside of the plasma. These 
assumptions are similar to those made in the statistical model [51] and take into account 
the in-medium modification of the particle states. 

In order then to calculate the phenomenologically relevant density of charmonium 
states at freeze out from the spectral function, we use the following tactic: We estimate 
hrst the contribution from the different states to the in-medium dilepton emission rate. 
From it we can compute the ratio of dileptons produced by T' and J/T. The ratio of 
the number densities of T' and J/T follows when we correct for the probability of the 
corresponding vacuum states to decay electromagnetically. 
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In detail, the dilepton emission rate (1.8) reads 


Ril oc 


dpQd?-p 


pjP) 

p2 


nsiPo)- 


(4.1) 


Here we use the fact that to leading order p depends only on = Vo ~ ^ind 
leave a more detailed analysis of the momentum dependence for future studies. After 
performing the change of variable oj = \/— p^, we obtain 

P^oc [ dud^p^^UBi^/u}^+ p‘^) . ^ (4.2) 

J + p2 

In this formula, the contribution from the different bound states will come from the 
corresponding peak area in p(a;)/a;^, see fig. 9. Hence we £t p(a;)/a;^ with a skewed 
Breit-Wigner (3.8) to distinguish the contribution from the different states n. From 
that fit we obtain the position and width of the spectral feature, i.e. the thermal mass 
Mn of the particle and its width. 

To calculate the integral (4.2) we checked numerically that it is possible to ap¬ 
proximate the Breit-Wigner peak by a delta peak, keeping the area A under the curve 
constant. Performing this, the integral (4.2) becomes 


Rjf (X A 


du d^pusi^/M^ + p2) 






(4.3) 


The fit of the charmonium spectrum in fig. 9 then yields 


R^- 

" = 0.023 ± 0.004. 


R 


.//'I' 

il 


(4.4) 


To obtain the number density we divide, as discussed, by the electromagnetic decay 
rate of a vacuum state, which is proportional to the square of the wave function 4) at 
r = 0 divided by the square of the mass of the state [52] 




N 


j/'i' 


T=Tc 




= 0.052 ± 0.009, 


(4.5) 


where we used the value at the origin of the zero temperature wave functions. A 
comparison of our result to those of the statistical model is shown in fig. 10. 

While the experimental determination of the T' to J/T ratio is challenging and 
currently limited by statistics, both ALICE [53] and CMS [54] have put forward first 


measurements of the double ratio 
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Figure 10. Comparison of predicted yields for and J/'I' from this work (orange) compared 
to predictions from the statistical model (green) and several experimental measurements. The 
plot structure and experimental data has been adapted from ref. [51]. 


systematic uncertainties makes this quantity the most robust measure of the in-medium 
relation between \h' and J/\h to date. Integrated over centrality CMS found that at 
the lowest accessible rapidities the ratio lies at 0.45 ± 0.13(stat) ± 0.07(syst) , while at 
larger rapidity 1.6 < \y\ < 2.4 it grows to values larger than one. 

Here we give a rough estimate of the double ratio by using experimental data ob¬ 
tained for prompt charmonium at y/s = 7TeV by the CMS collaboration. The rapidity 
averaged cross-section ratios for 4/' to J/4/ denoted by R in [55] are furthermore aver¬ 
aged over transverse momentum {R)pj. = 0.0378 ± 0.006, as the observed dependence 
on pj- is rather weak. To extract the number density ratio, the difference in branching 
fraction into dimuons must also be corrected, after which we obtain 




^=7TeV 

This leads to a double ratio of 

Nq/i 


, , BR(J/'^ u+u” , 

= (^)pT = 0.09 ± 0.015, 




(4.6) 


N 


j/'f 


N\i)i 


PbPb 


= 0.58 ±0.14 


(4.7) 


pp 


in which the errors in the individual contributions have been naively propagated under 
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Figure 11. charmonium (left) bottomonium (right) spectra at T = 3Tc from an extrapola¬ 
tion of the HTL fitted values for rriD- To visualize the uncertainties of the computation three 
curves according to the values of rriD and mo ^ 5m£) are shown. 


the assumption of being independent. Compared to the low rapidity measurements by 
CMS [54] we hnd good agreement within our still sizable error bars. 

4.2 Bottomonium at maximum LHC temperatures 

Given the small amount of bottom quarks produced in heavy-ion collisions, both at 
RHIC and LHC, one would naively expect a negligible amount of recombination, a 
thought which might however need reconsideration [56, 57]. The bottomonium states 
seen in lead-lead collisions are hence expected to be the ones that survive over the 
whole plasma evolution. A way to approximate their number is to read off the peak 
area from the spectrum at the highest temperature reached by the plasma. At LHC, 
where most experimental results for bottomonium have been obtained, the highest 
temperature according to ref. [58] corresponds to T ~ 3Tc from an analysis of harmonic 
yields. We hence consider the spectrum at T^ax = 3Tc in fig. 11. There, we see that 
only the ground state survives and is noticeably washed out. Of course not all the 
bottomonium formed will experience the highest temperature, one e.g. expects states 
to be also produced at the edge of the plasma, where the temperature is lower. 

If we assume the fate of the heavy quarkonium is determined at the maximum tem¬ 
perature reached by the plasma, we can estimate the ratio of bottomonium produced 
in lead-lead collisions to the amount of Bottom produced in proton-proton collisions, 
which is given by the ratio of the spectral weights at T^ax = 3Tc vs T = 0, 

^PbPb/pp = = 0.91 ± 0.02. (4.8) 

i r=oGT=o 

Here T denotes the width and C the amplitude of the corresponding Breit-Wigner type 
spectral feature. Note that the error bars might be underestimated as the Debye mass 
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had to be extrapolated far above the available temperatures using the HTL fit. The 
corresponding experimental results from CMS [59] are lower with i?pbPb/pp = 0.6 ±0.2. 

4.3 Melting temperatures 

As we found in fig. 4, the bound states disappear progressively with increasing T and 
we can attempt to define their melting point. In the case of a real potential, the 
melting temperature is straightforward to dehne from the disappearance of the purely 
real binding energy. For a complex potential the situation is more subtle, as the bound 
state broadens before it disappears. A popular choice is to dehne the melting of a state 
the moment its width equals its binding energy [10, 60]. 

Here we scan the spectrum at different temperatures, using the HTL based inter- 
and extrapolation of the values of the continuum corrected Debye mass niD- The 
bound state peak features are determined from a ht based on the skewed Lorentzian of 
eq.(3.8). The melting temperatures obtained in this way are given in tab. 6 and exhibit 
a hierarchical pattern, where only the ground states survive well above Tc^. 

Our observations are in qualitative agreement with earlier lattice QCD studies 
of charmonium correlators and spectral functions [62-70]. Note that most of these 
studies were performed in the quenched approximation or at rather large values of the 
light quark masses and so far no continuum extrapolation was performed. Results 
of spatial correlation functions and the corresponding spatial screening masses from 
lattice calculations in the charmoniun sector that include dynamical light quarks further 
support our hndings [71, 72]. Relativistic charmonium and bottomonium correlation 
functions computed in quenched QCD [70] have also shown that in the bottomonium 
channel thermal modifications set in at larger temperatures, well in the QGP phase, 
compared to the charmonium channel. The recent determination of bottomonium 
spectra in full QCD based on non-relativistic QCD [73, 74] corroborates an T(IS') 
ground state survival deep into the QGP. A thorough quantitative comparison of our 
lattice potential based spectra to direct lattice QGD computations will be part of a 
future study and is work in progress. 

A more precise determination of the T(IS') melting temperatures will require lattice 
simulations at higher temperatures, as for now we can only use an extrapolation of the 
HTL £t for mo, which becomes increasingly unreliable at high temperatures beyond 
T > 1.66Tc. 


®For charmonium a similar behavior has been anticipated in the in the sequential melting scenario 
of ref. [61]. 
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states 

J/T(15) 

T'(25) 

T(15) 

T(25) 

T(3F) 

T(4^) 

—i^’bind IT' 
-'melt I^C 

1 Q7+0.08 

-‘-•'J ' -0.07 

< 0.95 

2 

-‘-•^"^-0.05 

1 Q1+0.03 

< 0.95 


Table 6. Melting temperatures Tmeit of the different bound states defined by the point, at 
which the width of the state equates its binding energy calculated from Re[V] alone. The error 
on the determination of Tmeit takes into account the possible variation ofrriD as shown in fig. 3. 
Note that because of the lack of sufficient number of lattice ensembles below T = O.ObTc and 
the breakdown of the generalized Gauss law ansatz, we only give upper limits in this regime. 

5 Conclusion 

Heavy quarkonium is a vital probe to uncover the physics of the quark-gluon plasma 
created in heavy-ion collisions. Their non-relativistic nature opens up the possibility 
to describe their in-medium behavior by an effective potential entering a Schrodinger 
equation. We reported here the hrst quantitative phenomenological investigation of 
bottomonium and charmonium S-wave spectra at hnite T, based on the proper complex 
in-medium static potential extracted from Nf = 2 + 1 lattice QCD. 

To make possible the use of the lattice potential in a phenomenological application 
we deployed the generalized Gauss law ansatz, which provides analytic expressions 
for Re[V] and Im[V] that depend on a single temperature dependent parameter 
the Debye mass. Tuning rriD it was possible to reproduce the lattice values of Re[V] 
and to validate the up to now tentative values of Im[V] at all separation distances 
and temperatures investigated. After correcting for lattice discretization artifacts in 
the Debye mass we htted its temperature dependence successfully with a HTL based 
expression. 

To compute spectral functions for phenomenological inspection, we determined 
the vacuum parameters of the static potential in the continuum from a £t to the 
experimentally known S-, P- and D-wave bottomonium states. Finite T effects are 
incorporated though the continuum corrected Debye mass, extracted from the real-part 
of the in-medium lattice potential, leading to correctly renormalized hnite temperature 
expressions for Re[V] and Im[V] 

The central hndings in our adiabatic setup are a clear pattern of sequential melting 
of both bottomonium and charmonium with respect to their vacuum binding energies. 
The interplay between Debye screening and scattering with medium partons leads to 
characteristic shifts of the in-medium bound states to lower masses before they dissolve 
into the continuum. This effect is opposite to the relatively small gain in thermal mass. 

The availability of in-medium spectra with physical widths allowed us to estimate 
the T' to J/T ratio at the crossover transition temperature, relevant in the context 
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of the proposed statistical hadronization scenario. Our approach to compute the ratio 
of in-medium dimuon emission corrected by the emission rate of the corresponding 
vacuum states yields a value slightly larger than that of the statistical model but still 
much smaller than the experimental data m. p + p collisions. A rough estimate of the 
suppression of the T(IS') state in this adiabatic scenario on the other hand gave a value 
50% larger than experimentally observed. 

The static potential is only the starting point for a quantitatively reliable investiga¬ 
tion of heavy quarkonium. In particular for the lighter flavor charmonium the question 
of hnite mass effects in the potential should be addressed and efforts should be directed 
at determining the hrst correction Vi{r)/mQ at hnite temperature from the lattice. 
On the other hand direct lattice QCD studies, be it in the relativistic formulation or 
using the effective held theory NRQCD including velocity correction, will be essential 
to crosscheck the validity of the potential description. Comparisons with the data here 
can be made both on the level of spectra, as well as the corresponding Euclidean cor¬ 
relators. The necessity to solve an inherently ill-dehned problem in order to extract 
the spectral functions from lattice QCD current-current correlators directly remains a 
signihcant challenge. In order to surmount it both methodological progress in Bayesian 
and non-Bayesian spectral reconstruction approaches, as well as computational ehorts 
to generate lattices with more Euclidean time steps will be required. 

The lattice in-medium potential can also be used to setup a description of the 
real-time evolution of the heavy quarkonium wave function at hnite temperature in the 
context of open quantum systems. An approach based on a stochastic potential [19] 
appears promising, where a purely real potential, given by Re[V] is randomly perturbed 
at each time step by noise of strength related to Im[V]. If developed further it promises 
to provide vital and phenomenologically relevant information on e.g. the density matrix 
of states for the quarkonium system, which goes beyond what is accessible from within 
the spectral functions computed here. 

Directions for future work include the determination of the complex in-medium 
potential on lattices at higher temperature in order to avoid the extrapolations beyond 
T = 1.66T(7 required e.g. in the determination the T(IS') melting temperature in 
sec. 4.2. Another aspect of interest is the momentum dependence of the spectra which 
is also an active area of research in direct lattice QCD studies. 
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